#include <boost/range/concepts.hpp>
#include <iostream>
#include <vector>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/vector.h>

using DealiiSpMat = dealii::SparseMatrix<double>;
using DealiiVectorXd = dealii::Vector<double>;
using dealii::SparseDirectUMFPACK;
using dealii::SparsityPattern;
using std::cout;
using std::endl;
using u32 = unsigned int;
using std::vector;
void test_copy_matrix(DealiiSpMat &_Mat)
{
  cout << "test copy matrix in dealii:\n";
  const SparsityPattern& spM = _Mat.get_sparsity_pattern();
  DealiiSpMat A;
  SparsityPattern spA;
  spA.copy_from(spM);
  cout << "SparsityPattern spA is:\n";
  spA.print(std::cout);
  A.reinit(spA);
  A.copy_from(_Mat);
  // for(u32 i = 0; i < _Mat.m(); i++) {
  //   for(auto it = _Mat.begin(i); it != _Mat.end(i); it++) {
  //     A.add(i,it->column(),it->value());
  //   }
  // }
  cout << "SparseMatrix A is :\n";
  A.print(std::cout);
  cout << "Pattern of SparsityPattern A is :\n";
  A.print_pattern(std::cout);
}


void test_vec_push(DealiiSpMat& _Mat)
{
  vector<SparsityPattern> vec_spm;
  vector<DealiiSpMat> vec_m;
  vec_spm.reserve(3);
  vec_m.reserve(3);
  for(int i = 0 ; i < 3; i++){
    DealiiSpMat mm;
    SparsityPattern spmm;
    vec_m.push_back(mm);
    const SparsityPattern& spM = _Mat.get_sparsity_pattern();
    vec_spm.push_back(spmm);
    vec_spm.back().copy_from(spM);
    vec_m.back().reinit(vec_spm.back());
    for(u32 i = 0; i < _Mat.m(); i++) {
      for(auto it = _Mat.begin(i); it != _Mat.end(i); it++) {
        vec_m.back().add(i,it->column(),it->value());
      }
    }
  }
}


  int main()
  {
    {
      DealiiSpMat M;
      SparsityPattern SPM;
      SPM.reinit(5,5,5);
      for(u32 i = 0; i < 5; i++) {
        for(u32 j = i; j < 5; j++) {
          SPM.add(i,j);
        }
      }
      SPM.compress();
      M.reinit(SPM);
      cout << "SparsityPattern SPM is \n";
      SPM.print(std::cout);
      for(u32 i = 0; i < M.m(); i++) {
        for(auto it = M.begin(i); it != M.end(i); it++) {
          M.add(i,it->column(),1);
        }
      }
      cout << "SparseMatrix M is :\n";
      M.print(std::cout);
      test_copy_matrix(M);
    }
    {
      DealiiSpMat M;
      SparsityPattern SPM;
      SPM.reinit(5,5,5);
      for(u32 i = 0; i < 5; i++) {
        for(u32 j = i; j < 5; j++) {
          SPM.add(i,j);
        }
      }
      SPM.compress();
      M.reinit(SPM);
      for(u32 i = 0; i < M.m(); i++) {
        for(auto it = M.begin(i); it != M.end(i); it++) {
          M.add(i,it->column(),1);
        }
      }
      test_vec_push(M);
    }
  };
